home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 11 / Cream of the Crop 11-2.iso / os2 / rsynth1.zip / nsynth.c < prev    next >
Text File  |  1995-12-19  |  32KB  |  1,067 lines
  1. #include <config.h>
  2.  
  3.  
  4. /* $Id: nsynth.c,v 1.13 1994/11/08 13:30:50 a904209 Exp a904209 $
  5.  */
  6. char *nsynth_id = "$Id: nsynth.c,v 1.13 1994/11/08 13:30:50 a904209 Exp a904209 $";
  7.  
  8. /* Copyright            1982                    by Dennis H. Klatt
  9.  
  10.  *      Klatt synthesizer
  11.  *         Modified version of synthesizer described in
  12.  *         J. Acoust. Soc. Am., Mar. 1980. -- new voicing
  13.  *         source.
  14.  *
  15.  * Edit history
  16.  * 000001 10-Mar-83 DK  Initial creation.
  17.  * 000002  5-May-83 DK  Fix bug in operation of parallel F1
  18.  * 000003  7-Jul-83 DK  Allow parallel B1 to vary, and if ALL_PARALLEL,
  19.  *                      also allow B2 and B3 to vary
  20.  * 000004 26-Jul-83 DK  Get rid of mulsh, use short for VAX
  21.  * 000005 24-Oct-83 DK  Split off parwavtab.c, change short to int
  22.  * 000006 16-Nov-83 DK  Make samrate a variable, use exp(), cos() rand()
  23.  * 000007 17-Nov-83 DK  Convert to float, remove  cpsw, add set outsl
  24.  * 000008 28-Nov-83 DK  Add simple impulsive glottal source option
  25.  * 000009  7-Dec-83 DK  Use spkrdef[7] to select impulse or natural voicing
  26.  *                       and update cascade F1,..,F6 at update times
  27.  * 000010 19-Dec-83 DK  Add subroutine no_rad_char() to get rid of rad char
  28.  * 000011 28-Jan-84 DK  Allow up to 8 formants in cascade branch F7 fixed
  29.  *                       at 6.5 kHz, F8 fixed at 7.5 kHz
  30.  * 000012 14-Feb-84 DK  Fix bug in 'os' options so os>12 works
  31.  * 000013 17-May-84 DK  Add G0 code
  32.  * 000014 12-Mar-85 DHW modify for Haskins environment
  33.  * 000015 11-Jul-87 LG  modificiations for PC
  34.  * 000016 20-Apr-91 ATS Modified for SPARCSTATION
  35.  */
  36.  
  37. #include <useconfig.h>
  38. #include <stdio.h>
  39. #include <math.h>
  40. #include "proto.h"
  41. #include "nsynth.h"
  42. #ifndef PI
  43. #ifndef M_PI                      /* <math.h> */
  44. #define PI               3.1415927
  45. #else /* M_PI */
  46. #define PI               M_PI
  47. #endif /* M_PI */
  48. #endif
  49.  
  50. #ifdef __STDC__
  51. #define ONE 1.0F
  52. #else
  53. #define ONE 1.0
  54. #endif
  55.  
  56. typedef struct
  57.  {
  58.   char *name;
  59.   float a;
  60.   float b;
  61.   float c;
  62.   float p1;
  63.   float p2;
  64.  }
  65. resonator_t, *resonator_ptr;
  66.  
  67. /* Various global variables */
  68.  
  69. int time_count = 0;
  70. static warnsw;                    /* JPI added for f0 flutter */
  71.  
  72. /* COUNTERS */
  73.  
  74. static long nper;                 /* Current loc in voicing period   40000 samp/s */
  75.  
  76. /* COUNTER LIMITS */
  77.  
  78. static long T0;                   /* Fundamental period in output samples times 4 */
  79. static long nopen;                /* Number of samples in open phase of period  */
  80. static long nmod;                 /* Position in period to begin noise amp. modul */
  81.  
  82. /* Variables that have to stick around for awhile, and thus locals
  83.    are not appropriate 
  84.  */
  85.  
  86. /* Incoming parameter Variables which need to be updated synchronously  */
  87.  
  88. static long F0hz10;               /* Voicing fund freq in Hz  */
  89. static long F0hz10org;            /* Original Fundamental Freq */
  90. static long AVdb;                 /* Amp of voicing in dB,    0 to   70  */
  91. static long Kskew;                /* Skewness of alternate periods,0 to   40  */
  92.  
  93. /* Various amplitude variables used in main loop */
  94.  
  95. static float amp_voice;           /* AVdb converted to linear gain  */
  96. static float amp_bypas;           /* AB converted to linear gain  */
  97. static float par_amp_voice;       /* AVpdb converted to linear gain  */
  98. static float amp_aspir;           /* AP converted to linear gain  */
  99. static float amp_frica;           /* AF converted to linear gain  */
  100. static float amp_breth;           /* ATURB converted to linear gain  */
  101.  
  102. /* State variables of sound sources */
  103.  
  104. static long skew;                 /* Alternating jitter, in half-period units  */
  105.  
  106. static float natglot_a;           /* Makes waveshape of glottal pulse when open  */
  107. static float natglot_b;           /* Makes waveshape of glottal pulse when open  */
  108. static float vwave;               /* Ditto, but before multiplication by AVdb  */
  109. static float vlast;               /* Previous output of voice  */
  110. static float nlast;               /* Previous output of random number generator  */
  111. static float glotlast;            /* Previous value of glotout  */
  112. static float decay;               /* TLTdb converted to exponential time const  */
  113. static float onemd;               /* in voicing one-pole low-pass filter  */
  114. static float minus_pi_t;          /* func. of sample rate */
  115. static float two_pi_t;            /* func. of sample rate */
  116.  
  117.  
  118. /* INTERNAL MEMORY FOR DIGITAL RESONATORS AND ANTIRESONATOR  */
  119.  
  120. static resonator_t rnpp =
  121. {"parallel nasal pole"};
  122. static resonator_t r1p =
  123. {"parallel 1st formant"};
  124. static resonator_t r2p =
  125. {"parallel 2nd formant"};
  126. static resonator_t r3p =
  127. {"parallel 3rd formant"};
  128. static resonator_t r4p =
  129. {"parallel 4th formant"};
  130. static resonator_t r5p =
  131. {"parallel 5th formant"};
  132. static resonator_t r6p =
  133. {"parallel 6th formant"};
  134. static resonator_t r1c =
  135. {"cascade 1st formant"};
  136. static resonator_t r2c =
  137. {"cascade 2nd formant"};
  138. static resonator_t r3c =
  139. {"cascade 3rd formant"};
  140. static resonator_t r4c =
  141. {"cascade 4th formant"};
  142. static resonator_t r5c =
  143. {"cascade 5th formant"};
  144. static resonator_t r6c =
  145. {"cascade 6th formant"};
  146. static resonator_t r7c =
  147. {"cascade 7th formant"};
  148. static resonator_t r8c =
  149. {"cascade 8th formant"};
  150. static resonator_t rnpc =
  151. {"cascade nasal pole"};
  152. static resonator_t rnz =
  153. {"cascade nasal zero"};
  154. static resonator_t rgl =
  155. {"crit-damped glot low-pass filter"};
  156. static resonator_t rlp =
  157. {"downsamp low-pass filter"};
  158. static resonator_t rout =
  159. {"output low-pass"};
  160.  
  161. /*
  162.  * Constant natglot[] controls shape of glottal pulse as a function
  163.  * of desired duration of open phase N0
  164.  * (Note that N0 is specified in terms of 40,000 samples/sec of speech)
  165.  *
  166.  *    Assume voicing waveform V(t) has form: k1 t**2 - k2 t**3
  167.  *
  168.  *    If the radiation characterivative, a temporal derivative
  169.  *      is folded in, and we go from continuous time to discrete
  170.  *      integers n:  dV/dt = vwave[n]
  171.  *                         = sum over i=1,2,...,n of { a - (i * b) }
  172.  *                         = a n  -  b/2 n**2
  173.  *
  174.  *      where the  constants a and b control the detailed shape
  175.  *      and amplitude of the voicing waveform over the open
  176.  *      potion of the voicing cycle "nopen".
  177.  *
  178.  *    Let integral of dV/dt have no net dc flow --> a = (b * nopen) / 3
  179.  *
  180.  *    Let maximum of dUg(n)/dn be constant --> b = gain / (nopen * nopen)
  181.  *      meaning as nopen gets bigger, V has bigger peak proportional to n
  182.  *
  183.  *    Thus, to generate the table below for 40 <= nopen <= 263:
  184.  *
  185.  *      natglot[nopen - 40] = 1920000 / (nopen * nopen)
  186.  */
  187. static const short natglot[224] =
  188. {
  189.  1200, 1142, 1088, 1038, 991, 948, 907, 869, 833, 799,
  190.  768, 738, 710, 683, 658, 634, 612, 590, 570, 551,
  191.  533, 515, 499, 483, 468, 454, 440, 427, 415, 403,
  192.  391, 380, 370, 360, 350, 341, 332, 323, 315, 307,
  193.  300, 292, 285, 278, 272, 265, 259, 253, 247, 242,
  194.  237, 231, 226, 221, 217, 212, 208, 204, 199, 195,
  195.  192, 188, 184, 180, 177, 174, 170, 167, 164, 161,
  196.  158, 155, 153, 150, 147, 145, 142, 140, 137, 135,
  197.  133, 131, 128, 126, 124, 122, 120, 119, 117, 115,
  198.  113, 111, 110, 108, 106, 105, 103, 102, 100, 99,
  199.  97, 96, 95, 93, 92, 91, 90, 88, 87, 86,
  200.  85, 84, 83, 82, 80, 79, 78, 77, 76, 75,
  201.  75, 74, 73, 72, 71, 70, 69, 68, 68, 67,
  202.  66, 65, 64, 64, 63, 62, 61, 61, 60, 59,
  203.  59, 58, 57, 57, 56, 56, 55, 55, 54, 54,
  204.  53, 53, 52, 52, 51, 51, 50, 50, 49, 49,
  205.  48, 48, 47, 47, 46, 46, 45, 45, 44, 44,
  206.  43, 43, 42, 42, 41, 41, 41, 41, 40, 40,
  207.  39, 39, 38, 38, 38, 38, 37, 37, 36, 36,
  208.  36, 36, 35, 35, 35, 35, 34, 34, 33, 33,
  209.  33, 33, 32, 32, 32, 32, 31, 31, 31, 31,
  210.  30, 30, 30, 30, 29, 29, 29, 29, 28, 28,
  211.  28, 28, 27, 27
  212. };
  213.  
  214. /*
  215.  * Convertion table, db to linear, 87 dB --> 32767
  216.  *                                 86 dB --> 29491 (1 dB down = 0.5**1/6)
  217.  *                                 ...
  218.  *                                 81 dB --> 16384 (6 dB down = 0.5)
  219.  *                                 ...
  220.  *                                  0 dB -->     0
  221.  *
  222.  * The just noticeable difference for a change in intensity of a vowel
  223.  *   is approximately 1 dB.  Thus all amplitudes are quantized to 1 dB
  224.  *   steps.
  225.  */
  226.  
  227. static const float amptable[88] =
  228. {
  229.  0.0, 0.0, 0.0, 0.0, 0.0,
  230.  0.0, 0.0, 0.0, 0.0, 0.0,
  231.  0.0, 0.0, 0.0, 6.0, 7.0,
  232.  8.0, 9.0, 10.0, 11.0, 13.0,
  233.  14.0, 16.0, 18.0, 20.0, 22.0,
  234.  25.0, 28.0, 32.0, 35.0, 40.0,
  235.  45.0, 51.0, 57.0, 64.0, 71.0,
  236.  80.0, 90.0, 101.0, 114.0, 128.0,
  237.  142.0, 159.0, 179.0, 202.0, 227.0,
  238.  256.0, 284.0, 318.0, 359.0, 405.0,
  239.  455.0, 512.0, 568.0, 638.0, 719.0,
  240.  811.0, 911.0, 1024.0, 1137.0, 1276.0,
  241.  1438.0, 1622.0, 1823.0, 2048.0, 2273.0,
  242.  2552.0, 2875.0, 3244.0, 3645.0, 4096.0,
  243.  4547.0, 5104.0, 5751.0, 6488.0, 7291.0,
  244.  8192.0, 9093.0, 10207.0, 11502.0, 12976.0,
  245.  14582.0, 16384.0, 18350.0, 20644.0, 23429.0,
  246.  26214.0, 29491.0, 32767.0
  247. };
  248.  
  249. const char *par_name[] =
  250. {
  251.  "F0hz10",
  252.  "AVdb",
  253.  "F1hz", "B1hz",
  254.  "F2hz", "B2hz",
  255.  "F3hz", "B3hz",
  256.  "F4hz", "B4hz",
  257.  "F5hz", "B5hz",
  258.  "F6hz", "B6hz",
  259.  "FNZhz", "BNZhz",
  260.  "FNPhz", "BNPhz",
  261.  "AP",
  262.  "Kopen",
  263.  "Aturb",
  264.  "TLTdb",
  265.  "AF",
  266.  "Kskew",
  267.  "A1", "B1phz",
  268.  "A2", "B2phz",
  269.  "A3", "B3phz",
  270.  "A4", "B4phz",
  271.  "A5", "B5phz",
  272.  "A6", "B6phz",
  273.  "ANP",
  274.  "AB",
  275.  "AVpdb",
  276.  "Gain0"
  277. };
  278.  
  279. static void flutter PROTO((klatt_global_ptr globals, klatt_frame_ptr pars));
  280. static float resonator PROTO((resonator_ptr r, Float input));
  281. static float antiresonator PROTO((resonator_ptr r, Float input));
  282. static float impulsive_source PROTO((long nper));
  283. static float natural_source PROTO((long nper));
  284. static void setabc PROTO((long int f, long int bw, resonator_ptr rp));
  285. static void setabcg PROTO((long int f, long int bw, resonator_ptr rp, Float gain));
  286. static void setzeroabc PROTO((long int f, long int bw, resonator_ptr rp));
  287. static float DBtoLIN PROTO((klatt_global_ptr globals, long int dB));
  288. static float dBconvert PROTO((long int arg));
  289. static void overload_warning PROTO((klatt_global_ptr globals, long int arg));
  290. static short clip PROTO((klatt_global_ptr globals, Float input));
  291. static void pitch_synch_par_reset PROTO((klatt_global_ptr globals,
  292.                                          klatt_frame_ptr frame, long ns));
  293. static void frame_init PROTO((klatt_global_ptr globals, klatt_frame_ptr frame));
  294. void show_parms PROTO((klatt_global_ptr globals, int *pars));
  295.  
  296.  
  297. void
  298. show_parms(globals, pars)
  299. klatt_global_ptr globals;
  300. int *pars;
  301. {
  302.  int i;
  303.  static int names;
  304.  if ((names++ % 64) == 0)
  305.   {
  306.    for (i = 0; i < NPAR; i++)
  307.     printf("%s ", par_name[i]);
  308.    printf("\n");
  309.   }
  310.  for (i = 0; i < NPAR; i++)
  311.   {
  312.    printf("%*d ", (int) strlen(par_name[i]), pars[i]);
  313.   }
  314.  printf("\n");
  315. }
  316.  
  317. /*
  318.    function FLUTTER
  319.  
  320.    This function adds F0 flutter, as specified in:
  321.  
  322.    "Analysis, synthesis and perception of voice quality variations among
  323.    female and male talkers" D.H. Klatt and L.C. Klatt JASA 87(2) February 1990.
  324.    Flutter is added by applying a quasi-random element constructed from three
  325.    slowly varying sine waves.
  326.  */
  327. static void
  328. flutter(globals, pars)
  329. klatt_global_ptr globals;
  330. klatt_frame_ptr pars;
  331. {
  332.  long original_f0 = pars->F0hz10 / 10;
  333.  double fla = (double) globals->f0_flutter / 50;
  334.  double flb = (double) original_f0 / 100;
  335.  double flc = sin(2 * PI * 12.7 * time_count);
  336.  double fld = sin(2 * PI * 7.1 * time_count);
  337.  double fle = sin(2 * PI * 4.7 * time_count);
  338.  double delta_f0 = fla * flb * (flc + fld + fle) * 10;
  339.  F0hz10 += (long) delta_f0;
  340. }
  341.  
  342. static float
  343. impulsive_source(nper)
  344. long nper;
  345. {
  346.  static float doublet[] =
  347.  {0., 13000000., -13000000.};
  348.  if (nper < 3)
  349.   {
  350.    vwave = doublet[nper];
  351.   }
  352.  else
  353.   {
  354.    vwave = 0.0;
  355.   }
  356.  /* Low-pass filter the differenciated impulse with a critically-damped
  357.     second-order filter, time constant proportional to Kopen */
  358.  return resonator(&rgl, vwave);
  359. }
  360.  
  361.  
  362. /* Vwave is the differentiated glottal flow waveform, there is a weak
  363.    spectral zero around 800 Hz, magic constants a,b reset pitch-synch
  364.  */
  365.  
  366. static float
  367. natural_source(nper)
  368. long nper;
  369. {
  370.  float lgtemp;
  371.  /* See if glottis open */
  372.  if (nper < nopen)
  373.   {
  374.    natglot_a -= natglot_b;
  375.    vwave += natglot_a;
  376.    lgtemp = vwave * 0.028;        /* function of samp_rate ? */
  377.    return (lgtemp);
  378.   }
  379.  else
  380.   {
  381.    /* Glottis closed */
  382.    vwave = 0.0;
  383.    return (0.0);
  384.   }
  385. }
  386.  
  387. /*----------------------------------------------------------------------------*/
  388. /* Convert formant freqencies and bandwidth into
  389.    resonator difference equation coefficents
  390.  */
  391. static void
  392. setabc(f, bw, rp)
  393. long int f;                       /* Frequency of resonator in Hz  */
  394. long int bw;                      /* Bandwidth of resonator in Hz  */
  395. resonator_ptr rp;                 /* Are output coefficients  */
  396. {
  397.  double arg = minus_pi_t * bw;
  398.  float r = exp(arg);              /* Let r  =  exp(-pi bw t) */
  399.  rp->c = -(r * r);                /* Let c  =  -r**2 */
  400.  arg = two_pi_t * f;
  401.  rp->b = r * cos(arg) * 2.0;      /* Let b = r * 2*cos(2 pi f t) */
  402.  rp->a = 1.0 - rp->b - rp->c;     /* Let a = 1.0 - b - c */
  403. }
  404.  
  405. /* Convienience function for setting parallel resonators with gain */
  406. static void
  407. setabcg(f, bw, rp, gain)
  408. long int f;                       /* Frequency of resonator in Hz  */
  409. long int bw;                      /* Bandwidth of resonator in Hz  */
  410. resonator_ptr rp;                 /* Are output coefficients  */
  411. Float gain;
  412. {
  413.  setabc(f, bw, rp);
  414.  rp->a *= gain;
  415. }
  416.  
  417. /* Convert formant freqencies and bandwidth into
  418.  *      anti-resonator difference equation constants
  419.  */
  420. static void
  421. setzeroabc(f, bw, rp)
  422. long int f;                       /* Frequency of resonator in Hz  */
  423. long int bw;                      /* Bandwidth of resonator in Hz  */
  424. resonator_ptr rp;                 /* Are output coefficients  */
  425. {
  426.  setabc(f, bw, rp);               /* First compute ordinary resonator coefficients */
  427.  /* Now convert to antiresonator coefficients */
  428.  rp->a = 1.0 / rp->a;             /* a'=  1/a */
  429.  rp->b *= -rp->a;                 /* b'= -b/a */
  430.  rp->c *= -rp->a;                 /* c'= -c/a */
  431. }
  432.  
  433. /*----------------------------------------------------------------------------*/
  434.  
  435.  
  436. /* Convert from decibels to a linear scale factor */
  437. static float
  438. DBtoLIN(globals, dB)
  439. klatt_global_ptr globals;
  440. long int dB;
  441. {
  442.  /* Check limits or argument (can be removed in final product) */
  443.  if (dB < 0)
  444.   dB = 0;
  445.  else if (dB >= 88)
  446.   {
  447.    if (!globals->quiet_flag)
  448.     printf("Try to compute amptable[%ld]\n", dB);
  449.    dB = 87;
  450.   }
  451.  return amptable[dB] * 0.001;
  452. }
  453.  
  454. /* WHAT WERE THESE FOR ? */
  455. #define ACOEF           0.005
  456. #define BCOEF           (1.0 - ACOEF)    /* Slight decay to remove dc */
  457.  
  458. static float
  459. dBconvert(arg)
  460. long int arg;
  461. {
  462.  return 20.0 * log10((double) arg / 32767.0);
  463. }
  464.  
  465. static void
  466. overload_warning(globals, arg)
  467. klatt_global_ptr globals;
  468. long int arg;
  469. {
  470.  if (warnsw == 0)
  471.   {
  472.    warnsw++;
  473.    if (!globals->quiet_flag)
  474.     {
  475.      printf("\n* * * WARNING: ");
  476.      printf(" Signal at output of synthesizer (+%3.1f dB) exceeds 0 dB\n",
  477.             dBconvert(arg));
  478.     }
  479.   }
  480. }
  481.  
  482. /* Reset selected parameters pitch-synchronously */
  483.  
  484. static void
  485. pitch_synch_par_reset(globals, frame, ns)
  486. klatt_global_ptr globals;
  487. klatt_frame_ptr frame;
  488. long ns;
  489. {
  490.  long temp;
  491.  float temp1;
  492.  if (F0hz10 > 0)
  493.   {
  494.    T0 = (40 * globals->samrate) / F0hz10;
  495.    /* Period in samp*4 */
  496.    amp_voice = DBtoLIN(globals, AVdb);
  497.  
  498.    /* Duration of period before amplitude modulation */
  499.    nmod = T0;
  500.    if (AVdb > 0)
  501.     {
  502.      nmod >>= 1;
  503.     }
  504.  
  505.    /* Breathiness of voicing waveform */
  506.  
  507.    amp_breth = DBtoLIN(globals, frame->Aturb) * 0.1;
  508.  
  509.    /* Set open phase of glottal period */
  510.    /* where  40 <= open phase <= 263 */
  511.  
  512.    nopen = 4 * frame->Kopen;
  513.    if ((globals->glsource == IMPULSIVE) && (nopen > 263))
  514.     {
  515.      nopen = 263;
  516.     }
  517.  
  518.    if (nopen >= (T0 - 1))
  519.     {
  520.      nopen = T0 - 2;
  521.      if (!globals->quiet_flag)
  522.       {
  523.        printf("Warning: glottal open period cannot exceed T0, truncated\n");
  524.       }
  525.     }
  526.  
  527.    if (nopen < 40)
  528.     {
  529.      nopen = 40;                  /* F0 max = 1000 Hz */
  530.      if (!globals->quiet_flag)
  531.       {
  532.        printf("Warning: minimum glottal open period is 10 samples.\n");
  533.        printf("truncated, nopen = %ld\n", nopen);
  534.       }
  535.     }
  536.  
  537.    /* Reset a & b, which determine shape of "natural" glottal waveform */
  538.  
  539.    natglot_b = natglot[nopen - 40];
  540.    natglot_a = (natglot_b * nopen) * .333;
  541.  
  542.    /* Reset width of "impulsive" glottal pulse */
  543.  
  544.    temp = globals->samrate / nopen;
  545.    setabc(0L, temp, &rgl);
  546.  
  547.    /* Make gain at F1 about constant */
  548.  
  549.    temp1 = nopen * .00833;
  550.    rgl.a *= (temp1 * temp1);
  551.  
  552.    /* Truncate skewness so as not to exceed duration of closed phase
  553.       of glottal period */
  554.  
  555.    temp = T0 - nopen;
  556.    if (Kskew > temp)
  557.     {
  558.      if (!globals->quiet_flag)
  559.       {
  560.        printf("Kskew duration=%ld > glottal closed period=%ld, truncate\n",
  561.               Kskew, T0 - nopen);
  562.       }
  563.      Kskew = temp;
  564.     }
  565.    if (skew >= 0)
  566.     {
  567.      skew = Kskew;                /* Reset skew to requested Kskew */
  568.     }
  569.    else
  570.     {
  571.      skew = -Kskew;
  572.     }
  573.  
  574.    /* Add skewness to closed portion of voicing period */
  575.  
  576.    T0 = T0 + skew;
  577.    skew = -skew;
  578.   }
  579.  else
  580.   {
  581.    T0 = 4;                        /* Default for f0 undefined */
  582.    amp_voice = 0.0;
  583.    nmod = T0;
  584.    amp_breth = 0.0;
  585.    natglot_a = 0.0;
  586.    natglot_b = 0.0;
  587.   }
  588.  
  589.  /* Reset these pars pitch synchronously or at update rate if f0=0 */
  590.  
  591.  if ((T0 != 4) || (ns == 0))
  592.   {
  593.    /* Set one-pole low-pass filter that tilts glottal source */
  594.    decay = (0.033 * frame->TLTdb);    /* Function of samp_rate ? */
  595.    if (decay > 0.0)
  596.     {
  597.      onemd = 1.0 - decay;
  598.     }
  599.    else
  600.     {
  601.      onemd = 1.0;
  602.     }
  603.   }
  604. }
  605.  
  606. /* Get variable parameters from host computer,
  607.    initially also get definition of fixed pars
  608.  */
  609.  
  610. static void
  611. frame_init(globals, frame)
  612. klatt_global_ptr globals;
  613. klatt_frame_ptr frame;
  614. {
  615.  long Gain0;                      /* Overall gain, 60 dB is unity  0 to   60  */
  616.  float amp_parF1;                 /* A1 converted to linear gain  */
  617.  float amp_parFN;                 /* ANP converted to linear gain  */
  618.  float amp_parF2;                 /* A2 converted to linear gain  */
  619.  float amp_parF3;                 /* A3 converted to linear gain  */
  620.  float amp_parF4;                 /* A4 converted to linear gain  */
  621.  float amp_parF5;                 /* A5 converted to linear gain  */
  622.  float amp_parF6;                 /* A6 converted to linear gain  */
  623.  
  624. #if DEBUG
  625.  show_parms(globals, (int *) frame);
  626. #endif
  627.  
  628.  /*
  629.     Read  speech frame definition into temp store
  630.     and move some parameters into active use immediately
  631.     (voice-excited ones are updated pitch synchronously
  632.     to avoid waveform glitches).
  633.   */
  634.  
  635.  F0hz10 = frame->F0hz10;
  636.  F0hz10org = F0hz10; 
  637.  AVdb = frame->AVdb - 7;
  638.  if (AVdb < 0)
  639.   AVdb = 0;
  640.  
  641.  amp_aspir = DBtoLIN(globals, frame->ASP) * .05;
  642.  amp_frica = DBtoLIN(globals, frame->AF) * 0.25;
  643.  
  644.  Kskew = frame->Kskew;
  645.  par_amp_voice = DBtoLIN(globals, frame->AVpdb);
  646.  
  647.  /* Fudge factors (which comprehend affects of formants on each other?)
  648.     with these in place ALL_PARALLEL should sound as close as 
  649.     possible to CASCADE_PARALLEL.
  650.     Possible problem feeding in Holmes's amplitudes given this.
  651.   */
  652.  amp_parF1 = DBtoLIN(globals, frame->A1) * 0.4;    /* -7.96 dB */
  653.  amp_parF2 = DBtoLIN(globals, frame->A2) * 0.15;    /* -16.5 dB */
  654.  amp_parF3 = DBtoLIN(globals, frame->A3) * 0.06;    /* -24.4 dB */
  655.  amp_parF4 = DBtoLIN(globals, frame->A4) * 0.04;    /* -28.0 dB */
  656.  amp_parF5 = DBtoLIN(globals, frame->A5) * 0.022;    /* -33.2 dB */
  657.  amp_parF6 = DBtoLIN(globals, frame->A6) * 0.03;    /* -30.5 dB */
  658.  amp_parFN = DBtoLIN(globals, frame->ANP) * 0.6;    /* -4.44 dB */
  659.  amp_bypas = DBtoLIN(globals, frame->AB) * 0.05;    /* -26.0 db */
  660.  
  661.  if (globals->nfcascade >= 8)
  662.   {
  663.    /* Inside Nyquist rate ? */
  664.    if (globals->samrate >= 16000)
  665.     setabc(7500, 600, &r8c);
  666.    else
  667.     globals->nfcascade = 6;
  668.   }
  669.  
  670.  if (globals->nfcascade >= 7)
  671.   {
  672.    /* Inside Nyquist rate ? */
  673.    if (globals->samrate >= 16000)
  674.     setabc(6500, 500, &r7c);
  675.    else
  676.     globals->nfcascade = 6;
  677.   }
  678.  
  679.  /* Set coefficients of variable cascade resonators */
  680.  
  681.  if (globals->nfcascade >= 6)
  682.   setabc(frame->F6hz, frame->B6hz, &r6c);
  683.  
  684.  if (globals->nfcascade >= 5)
  685.   setabc(frame->F5hz, frame->B5hz, &r5c);
  686.  
  687.  setabc(frame->F4hz, frame->B4hz, &r4c);
  688.  setabc(frame->F3hz, frame->B3hz, &r3c);
  689.  setabc(frame->F2hz, frame->B2hz, &r2c);
  690.  setabc(frame->F1hz, frame->B1hz, &r1c);
  691.  
  692.  /* Set coeficients of nasal resonator and zero antiresonator */
  693.  setabc(frame->FNPhz, frame->BNPhz, &rnpc);
  694.  setzeroabc(frame->FNZhz, frame->BNZhz, &rnz);
  695.  
  696.  /* Set coefficients of parallel resonators, and amplitude of outputs */
  697.  setabcg(frame->F1hz, frame->B1phz, &r1p, amp_parF1);
  698.  setabcg(frame->FNPhz, frame->BNPhz, &rnpp, amp_parFN);
  699.  setabcg(frame->F2hz, frame->B2phz, &r2p, amp_parF2);
  700.  setabcg(frame->F3hz, frame->B3phz, &r3p, amp_parF3);
  701.  setabcg(frame->F4hz, frame->B4phz, &r4p, amp_parF4);
  702.  setabcg(frame->F5hz, frame->B5phz, &r5p, amp_parF5);
  703.  setabcg(frame->F6hz, frame->B6phz, &r6p, amp_parF6);
  704.  
  705.  
  706.  /* fold overall gain into output resonator */
  707.  Gain0 = frame->Gain0 - 3;
  708.  if (Gain0 <= 0)
  709.   Gain0 = 57;
  710.  /* output low-pass filter - resonator with freq 0 and BW = globals->samrate
  711.     Thus 3db point is globals->samrate/2 i.e. Nyquist limit.
  712.     Only 3db down seems rather mild...
  713.   */
  714.  setabcg(0L, (long) globals->samrate, &rout, DBtoLIN(globals, Gain0));
  715. }
  716.  
  717. static short
  718. clip(globals, input)
  719. klatt_global_ptr globals;
  720. Float input;
  721. {
  722.  long temp = input;
  723.  /* clip on boundaries of 16-bit word */
  724.  if (temp < -32767)
  725.   {
  726.    overload_warning(globals, -temp);
  727.    temp = -32767;
  728.   }
  729.  else if (temp > 32767)
  730.   {
  731.    overload_warning(globals, temp);
  732.    temp = 32767;
  733.   }
  734.  return (temp);
  735. }
  736.  
  737. /* Generic resonator function */
  738. static float
  739. resonator(r, input)
  740. resonator_ptr r;
  741. Float input;
  742. {
  743.  register float x ;
  744. /*
  745.  if(r->p1 < 1e-25 && r->p1 > -1e-25)
  746.     r->p1 = 0.0;
  747.  if(r->p2 < 1e-25 && r->p2 > -1e-25)
  748.     r->p2 = 0.0;
  749. */
  750.  x = r->a * input + r->b * r->p1 + r->c * r->p2;
  751.  r->p2 = r->p1;
  752.  r->p1 = x;
  753.  return x;
  754. }
  755.  
  756. /* Generic anti-resonator function
  757.    Same as resonator except that a,b,c need to be set with setzeroabc()
  758.    and we save inputs in p1/p2 rather than outputs.
  759.    There is currently only one of these - "rnz"
  760.  */
  761. /*  Output = (rnz.a * input) + (rnz.b * oldin1) + (rnz.c * oldin2) */
  762.  
  763. static float
  764. antiresonator(r, input)
  765. resonator_ptr r;
  766. Float input;
  767. {
  768.  register float x;
  769.  
  770. /*
  771.  if(r->p1 < 1e-25 && r->p1 > -1e-25)
  772.     r->p1 = 0.0;
  773.  if(r->p2 < 1e-25 && r->p2 > -1e-25)
  774.     r->p2 = 0.0;
  775. */
  776.  x = r->a * input + r->b * r->p1 + r->c * r->p2;
  777.  r->p2 = r->p1;
  778.  r->p1 = input;
  779.  return x;
  780. }
  781.  
  782. /*
  783.    function PARWAV
  784.  
  785.    CONVERT FRAME OF PARAMETER DATA TO A WAVEFORM CHUNK
  786.    Synthesize globals->nspfr samples of waveform and store in jwave[].
  787.  */
  788.  
  789. void
  790. parwave(globals, frame, jwave)
  791. klatt_global_ptr globals;
  792. klatt_frame_ptr frame;
  793. short int *jwave;
  794. {
  795.  long ns;
  796.  float out = 0.0;
  797.  /* Output of cascade branch, also final output  */
  798.  
  799.  /* Initialize synthesizer and get specification for current speech
  800.     frame from host microcomputer */
  801.  
  802.  frame_init(globals, frame);
  803.  
  804.  if (globals->f0_flutter != 0)
  805.   {
  806.    time_count++;                  /* used for f0 flutter */
  807.    flutter(globals, frame);       /* add f0 flutter */
  808.   }
  809.  
  810.  /* MAIN LOOP, for each output sample of current frame: */
  811.  
  812.  for (ns = 0; ns < globals->nspfr; ns++)
  813.   {
  814.    static unsigned long seed = 5; /* Fixed staring value */
  815.    float noise;
  816.    int n4;
  817.    float sourc;                   /* Sound source if all-parallel config used  */
  818.    float glotout;                 /* Output of glottal sound source  */
  819.    float par_glotout;             /* Output of parallelglottal sound sourc  */
  820.    float voice;                   /* Current sample of voicing waveform  */
  821.    float frics;                   /* Frication sound source  */
  822.    float aspiration;              /* Aspiration sound source  */
  823.    long nrand;                    /* Varible used by random number generator  */
  824.  
  825.    /* Our own code like rand(), but portable
  826.       whole upper 31 bits of seed random 
  827.       assumes 32-bit unsigned arithmetic
  828.       with untested code to handle larger.
  829.     */
  830.    seed = seed * 1664525 + 1;
  831.    if (8 * sizeof(unsigned long) > 32)
  832.          seed &= 0xFFFFFFFF;
  833.  
  834.    /* Shift top bits of seed up to top of long then back down to LS 14 bits */
  835.    /* Assumes 8 bits per sizeof unit i.e. a "byte" */
  836.    nrand = (((long) seed) << (8 * sizeof(long) - 32)) >> (8 * sizeof(long) - 14);
  837.  
  838.    /* Tilt down noise spectrum by soft low-pass filter having
  839.     *    a pole near the origin in the z-plane, i.e.
  840.     *    output = input + (0.75 * lastoutput) */
  841.  
  842.    noise = nrand + (0.75 * nlast);    /* Function of samp_rate ? */
  843.    nlast = noise;
  844.  
  845.    /* Amplitude modulate noise (reduce noise amplitude during
  846.       second half of glottal period) if voicing simultaneously present
  847.     */
  848.  
  849.    if (nper > nmod)
  850.     {
  851.      noise *= 0.5;
  852.     }
  853.  
  854.    /* Compute frication noise */
  855.    sourc = frics = amp_frica * noise;
  856.  
  857.    /* Compute voicing waveform : (run glottal source simulation at
  858.       4 times normal sample rate to minimize quantization noise in 
  859.       period of female voice)
  860.     */
  861.  
  862.    for (n4 = 0; n4 < 4; n4++)
  863.     {
  864.      if (globals->glsource == IMPULSIVE)
  865.       {
  866.        /* Use impulsive glottal source */
  867.        voice = impulsive_source(nper);
  868.       }
  869.      else
  870.       {
  871.        /* Or use a more-natural-shaped source waveform with excitation
  872.           occurring both upon opening and upon closure, stronest at closure */
  873.        voice = natural_source(nper);
  874.       }
  875.  
  876.      /* Reset period when counter 'nper' reaches T0 */
  877.      if (nper >= T0)
  878.       {
  879.        nper = 0;
  880.        pitch_synch_par_reset(globals, frame, ns);
  881.       }
  882.  
  883.      /* Low-pass filter voicing waveform before downsampling from 4*globals->samrate */
  884.      /* to globals->samrate samples/sec.  Resonator f=.09*globals->samrate, bw=.06*globals->samrate  */
  885.  
  886.      voice = resonator(&rlp, voice);    /* in=voice, out=voice */
  887.  
  888.      /* Increment counter that keeps track of 4*globals->samrate samples/sec */
  889.      nper++;
  890.     }
  891.  
  892.    /* Tilt spectrum of voicing source down by soft low-pass filtering, amount
  893.       of tilt determined by TLTdb
  894.     */
  895.    voice = (voice * onemd) + (vlast * decay);
  896.    vlast = voice;
  897.  
  898.    /* Add breathiness during glottal open phase */
  899.    if (nper < nopen)
  900.     {
  901.      /* Amount of breathiness determined by parameter Aturb */
  902.      /* Use nrand rather than noise because noise is low-passed */
  903.      voice += amp_breth * nrand;
  904.     }
  905.  
  906.    /* Set amplitude of voicing */
  907.  
  908.    glotout = amp_voice * voice;
  909.  
  910.    /* Compute aspiration amplitude and add to voicing source */
  911.    aspiration = amp_aspir * noise;
  912.    glotout += aspiration;
  913.  
  914.    par_glotout = glotout;
  915.  
  916.    if (globals->synthesis_model != ALL_PARALLEL)
  917.     {
  918.      /* Cascade vocal tract, excited by laryngeal sources.
  919.         Nasal antiresonator, then formants FNP, F5, F4, F3, F2, F1
  920.       */
  921.      float rnzout = antiresonator(&rnz, glotout);    /* Output of cascade nazal zero resonator  */
  922.      float casc_next_in = resonator(&rnpc, rnzout);    /* in=rnzout, out=rnpc.p1 */
  923.  
  924.      /* Recoded from sequence of if's to use C's fall through switch
  925.         semantics. May allow compiler to optimize
  926.       */
  927.      switch (globals->nfcascade)
  928.       {
  929.        case 8:
  930.         casc_next_in = resonator(&r8c, casc_next_in);    /* Do not use unless samrat = 16000 */
  931.        case 7:
  932.         casc_next_in = resonator(&r7c, casc_next_in);    /* Do not use unless samrat = 16000 */
  933.        case 6:
  934.         casc_next_in = resonator(&r6c, casc_next_in);    /* Do not use unless long vocal tract or samrat increased */
  935.        case 5:
  936.         casc_next_in = resonator(&r5c, casc_next_in);
  937.        case 4:
  938.         casc_next_in = resonator(&r4c, casc_next_in);
  939.        case 3:
  940.         casc_next_in = resonator(&r3c, casc_next_in);
  941.        case 2:
  942.         casc_next_in = resonator(&r2c, casc_next_in);
  943.        case 1:
  944.         out = resonator(&r1c, casc_next_in);
  945.         break;
  946.        default:
  947.         out = 0.0;
  948.       }
  949. #if 0
  950.      /* Excite parallel F1 and FNP by voicing waveform */
  951.      /* Source is voicing plus aspiration */
  952.      /* Add in phase, boost lows for nasalized */
  953.      out += (resonator(&rnpp, par_glotout) + resonator(&r1p, par_glotout));
  954. #endif
  955.     }
  956.    else
  957.     {
  958.      /* Is ALL_PARALLEL */
  959.      /* NIS - rsynth "hack"
  960.         As Holmes' scheme is weak at nasals and (physically) nasal cavity
  961.         is "back near glottis" feed glottal source through nasal resonators
  962.         Don't think this is quite right, but improves things a bit
  963.       */
  964.      par_glotout = antiresonator(&rnz, par_glotout);
  965.      par_glotout = resonator(&rnpc, par_glotout);
  966.      /* And just use r1p NOT rnpp */
  967.      out = resonator(&r1p, par_glotout);
  968.      /* Sound sourc for other parallel resonators is frication
  969.         plus first difference of voicing waveform.
  970.       */
  971.      sourc += (par_glotout - glotlast);
  972.      glotlast = par_glotout;
  973.     }
  974.  
  975.    /* Standard parallel vocal tract
  976.       Formants F6,F5,F4,F3,F2, outputs added with alternating sign
  977.     */
  978.    out = resonator(&r6p, sourc) - out;
  979.    out = resonator(&r5p, sourc) - out;
  980.    out = resonator(&r4p, sourc) - out;
  981.    out = resonator(&r3p, sourc) - out;
  982.    out = resonator(&r2p, sourc) - out;
  983.  
  984.    out = amp_bypas * sourc - out;
  985.  
  986.    out = resonator(&rout, out);
  987.    *jwave++ = clip(globals, out); /* Convert back to integer */
  988.   }
  989. }
  990.  
  991. void
  992. parwave_init(globals)
  993. klatt_global_ptr globals;
  994. {
  995.  long FLPhz = (950 * globals->samrate) / 10000;
  996.  long BLPhz = (630 * globals->samrate) / 10000;
  997.  
  998.  minus_pi_t = -PI / globals->samrate;
  999.  two_pi_t = -2.0 * minus_pi_t;
  1000.  
  1001.  setabc(FLPhz, BLPhz, &rlp);
  1002.  nper = 0;                        /* LG */
  1003.  T0 = 0;                          /* LG */
  1004.  
  1005.  rnpp.p1 = 0;                     /* parallel nasal pole  */
  1006.  rnpp.p2 = 0;
  1007.  
  1008.  r1p.p1 = 0;                      /* parallel 1st formant */
  1009.  r1p.p2 = 0;
  1010.  
  1011.  r2p.p1 = 0;                      /* parallel 2nd formant */
  1012.  r2p.p2 = 0;
  1013.  
  1014.  r3p.p1 = 0;                      /* parallel 3rd formant */
  1015.  r3p.p2 = 0;
  1016.  
  1017.  r4p.p1 = 0;                      /* parallel 4th formant */
  1018.  r4p.p2 = 0;
  1019.  
  1020.  r5p.p1 = 0;                      /* parallel 5th formant */
  1021.  r5p.p2 = 0;
  1022.  
  1023.  r6p.p1 = 0;                      /* parallel 6th formant */
  1024.  r6p.p2 = 0;
  1025.  
  1026.  r1c.p1 = 0;                      /* cascade 1st formant  */
  1027.  r1c.p2 = 0;
  1028.  
  1029.  r2c.p1 = 0;                      /* cascade 2nd formant  */
  1030.  r2c.p2 = 0;
  1031.  
  1032.  r3c.p1 = 0;                      /* cascade 3rd formant  */
  1033.  r3c.p2 = 0;
  1034.  
  1035.  r4c.p1 = 0;                      /* cascade 4th formant  */
  1036.  r4c.p2 = 0;
  1037.  
  1038.  r5c.p1 = 0;                      /* cascade 5th formant  */
  1039.  r5c.p2 = 0;
  1040.  
  1041.  r6c.p1 = 0;                      /* cascade 6th formant  */
  1042.  r6c.p2 = 0;
  1043.  
  1044.  r7c.p1 = 0;
  1045.  r7c.p2 = 0;
  1046.  
  1047.  r8c.p1 = 0;
  1048.  r8c.p2 = 0;
  1049.  
  1050.  rnpc.p1 = 0;                     /* cascade nasal pole  */
  1051.  rnpc.p2 = 0;
  1052.  
  1053.  rnz.p1 = 0;                      /* cascade nasal zero  */
  1054.  rnz.p2 = 0;
  1055.  
  1056.  rgl.p1 = 0;                      /* crit-damped glot low-pass filter */
  1057.  rgl.p2 = 0;
  1058.  
  1059.  rlp.p1 = 0;                      /* downsamp low-pass filter  */
  1060.  rlp.p2 = 0;
  1061.  
  1062.  vlast = 0;                       /* Previous output of voice  */
  1063.  nlast = 0;                       /* Previous output of random number generator  */
  1064.  glotlast = 0;                    /* Previous value of glotout  */
  1065.  warnsw = 0;
  1066. }
  1067.